1 Introduction

The Rise 2 Flow study used a number of instruments during the entry, exit, and weekly sessions (additionally personal data was collected during the entry session).

For entry and exit sessions six instruments were used:

  • Mindful attention awareness scale (MAAS).1
  • Scale of positive and negative experience (SPANE).2
  • Personal well-being index (PWB).3
  • Positive thinking (PST).4
  • Self efficacy (SE).5
  • Perceived productivity (PP).6

In addition, a seventh instrument, mini-IPIP, was added to the exit session to control for personality type.7

For the weekly sessions the WHO-5 well-being index was used (WHO-5).8

To summarize, the study was executed according to,

  1. An entry survey, at time \(t_0\) (MAAS, etc.)
  2. A weekly survey, at time \(t_{{t_0 + 1}}, \ldots, t_{{t_0+13}}\)
  3. Daily surveys
  4. An exit survey, at time \(t_{0+11}\) (MAAS, etc. and IPIP)

The main question the study tries to answer is: Is there a difference in the responses (i.e., higher/lower values on responses), across all or among some instruments, over time?

In order to answer that question we can first of all compare outcomes of \(t_0\) and \(t_1\), and contrast with personal data, e.g., age, gender, IPIP. Next, the weekly session can be used to investigate the trend over time, and perhaps see where a noticeable effect, if it exists, shows up.

2 Data cleaning

First we load all five datasets. For each dataset, make sure we set correct variable types.

entry <- read.xlsx(here("data/R2F Quant Comp Data Entry.xlsx")) # all instruments
exit <- read.xlsx(here("data/R2F Quant Comp Data Exit.xlsx")) # same here
weekly <- read.xlsx(here("data/R2F Quant Weekly.xlsx"), detectDates = T)
ipip <- read.xlsx(here("data/R2F IPIP.xlsx"))
pers <- read.xlsx(here("data/R2F Pers Data.xlsx")) # personal data
daily <- read.xlsx(here("data/R2F Daily.xlsx"), detectDates = T)

2.1 entry and exit datasets

We need to set a number of columns as factors (ordered factors in many cases). Additionally, when we have Yes/No answers, as we do in some questions below, we set them to 0/1 since they will be modeled with a Bernoulli likelihood.

First we clean the entry dataset.

for(j in 3:17){
  set(entry, i=NULL, j=j, value=factor(entry[[j]], levels = c("Almost never", "Very infrequently", "Somewhat infrequently", "Somewhat frequently", "Very frequently", "Almost always"), ordered = TRUE))
}

for(j in 18:29){
  set(entry, i=NULL, j=j, value=factor(entry[[j]], levels = c("Very rarely or never", "Rarely", "Sometimes", "Often", "Very often or always"), ordered = TRUE))
}

for(j in 30:37){
  set(entry, i=NULL, j=j, value=factor(entry[[j]], levels = c("Strongly disagree", "Disagree", "Slightly disagree", "Mixed or neither agree nor disagree", "Slightly agree", "Agree", "Strongly agree"), ordered = TRUE))
}

for(i in 38:59)
  entry[,eval(i)] <- ifelse(entry[eval(i)] == 'Yes', 1, 0)

for(j in 60:69){
  set(entry, i=NULL, j=j, value=factor(entry[[j]], levels = c("Hardly true","Rather true","Exactly true"), ordered = TRUE))
}

for(j in 70:76){
  set(entry, i=NULL, j=j, value=factor(entry[[j]], levels = c("None of the time","A little of the time","Some of the time","Most of the time","All of the time"), ordered = TRUE))
}

for(j in 77:79){
  set(entry, i=NULL, j=j, value=factor(entry[[j]], levels = c("1","2","3","4","5","6","7","8","9","10"), ordered = TRUE))
}

for(j in 80:80){
  set(entry, i=NULL, j=j, value=factor(entry[[j]], levels = c("You were a lot worse than other workers","You were a little worse than other workers","You were somewhat worse than other workers","You were about average","You were somewhat better than other workers","You were a little better than other workers","You were a lot better than other workers"), ordered = TRUE))
}

Next, we clean the exit dataset, which consists of the same variables as above.

for(j in 3:17){
  set(exit, i=NULL, j=j, value=factor(exit[[j]], levels = c("Almost never", "Very infrequently", "Somewhat infrequently", "Somewhat frequently", "Very frequently", "Almost always"), ordered = TRUE))
}

for(j in 18:29){
  set(exit, i=NULL, j=j, value=factor(exit[[j]], levels = c("Very rarely or never", "Rarely", "Sometimes", "Often", "Very often or always"), ordered = TRUE))
}

for(j in 30:37){
  set(exit, i=NULL, j=j, value=factor(exit[[j]], levels = c("Strongly disagree", "Disagree", "Slightly disagree", "Mixed or neither agree nor disagree", "Slightly agree", "Agree", "Strongly agree"), ordered = TRUE))
}

for(i in 38:59)
  exit[,eval(i)] <- ifelse(exit[eval(i)] == 'Yes', 1, 0)


for(j in 60:69){
  set(exit, i=NULL, j=j, value=factor(exit[[j]], levels = c("Hardly true","Rather true","Exactly true"), ordered = TRUE))
}

for(j in 70:76){
  set(exit, i=NULL, j=j, value=factor(exit[[j]], levels = c("None of the time","A little of the time","Some of the time","Most of the time","All of the time"), ordered = TRUE))
}

for(j in 77:79){
  set(exit, i=NULL, j=j, value=factor(exit[[j]], levels = c("1","2","3","4","5","6","7","8","9","10"), ordered = TRUE))
}

for(j in 80:80){
  set(exit, i=NULL, j=j, value=factor(exit[[j]], levels = c("You were a lot worse than other workers","You were a little worse than other workers","You were somewhat worse than other workers","You were about average","You were somewhat better than other workers","You were a little better than other workers","You were a lot better than other workers"), ordered = TRUE))
}

2.2 weekly dataset

The weekly dataset contains five columns of the same type, which we set as ordered categorical.

for(j in 4:8){
  set(weekly, i=NULL, j=j, value=factor(weekly[[j]], levels = c("At no time", "Less than half of the time", "Some of the time", "More than half of the time",  "Most of the time", "All of the time"), ordered = TRUE))
}

2.3 ipip dataset

The mini-IPIP instrument uses only five levels for all 20 questions (ordered categorical).

for(j in 3:22){
  set(ipip, i=NULL, j=j, value=factor(ipip[[j]], levels = c("Disagree", "Somewhat disagree", "Neither agree nor disagree", "Somewhat agree", "Agree"), ordered = TRUE))
}

2.4 pers dataset

In the pers dataset we have a mix of variable types we need to deal with. Age should be scaled, Gender and Occupation set to 1/2, and Living conditions set to factors (unordered).

# set 1/2 as indicators for gender
pers$Gender_0_1 <- ifelse(pers$Gender_0_1 == "Man/Transman", 1, 2)

# many categories exist for occupation, but let's simplify this so we can 
# use this variable for something. Indicate if they are students or not.
pers$Occupation <- ifelse(pers$Occupation == "Student", 1, 2)

# make living condition a categorical factor (unordered)
pers$Living_1_4 <- factor(pers$Living_1_4, levels = c("I live by myself", "I live in shared housing", "I live with a partner", "I live with my family"))

# remove the last three columns since we have >64% NAs
pers[,7:9] <- NULL

2.5 Missingness analysis

Check how many cells are NA.

table(is.na(entry))
## 
## FALSE  TRUE 
##  6862    98
table(is.na(exit))
## 
## FALSE  TRUE 
##  2684    36
table(is.na(weekly))
## 
## FALSE  TRUE 
##  2390     2
table(is.na(pers))
## 
## FALSE  TRUE 
##   519     3

Clearly we have missing data here, but it’s only a few percentages (\(0.6\)–\(1.8\)%). Let’s use complete case analysis for now (i.e., remove all rows containing NAs). Later, if needed, we can model the missingness in a principled Bayesian way.

entry <- entry[complete.cases(entry), ]
exit <- exit[complete.cases(exit), ]
weekly <- weekly[complete.cases(weekly), ]
pers <- pers[complete.cases(pers), ]

Remember, the main question the study tries to answer is: Is there a difference in the responses (i.e., higher/lower values on responses), across all or among some instruments, over time? This means that for now we can discard of the weekly data and instead focus on the entry and exit instruments, together with personal data (i.e., age, gender, occupation, and living conditions).

In order for us to conduct inferences along this line of thinking, we will make use of dummy variable regression estimators (DVRE), which is numerically identical to deviation score estimators. The DVRE approach dummy encodes our time variable \(t\) and sets an index \(0/1\), where \(t_0 = 0\) and \(t_1 = 1\). In short, each subject (ID) will have two rows where one row is the entry instrument \(t = 0\), and one row has the exit instrument \(t = 1\). We’ll then see if our predictors are useful for predicting the outcome (i.e., is there a significant difference in the \(\beta\) estimators) and then also investigate if the outcome at \(t_0\) equals \(t_1\) (and if they aren’t the same, how much do they differ?)

# convert Age to numeric and standardize (note we have two NAs, so we suppress the warning we get)
suppressWarnings(pers$Age <- scale(as.numeric(pers$Age)))

# create new column t, and set to 0 or 1
entry$t <- 0
exit$t <- 1

# Remove all suffixes in both files so we have the same column names
colnames(entry) <- gsub('_en', '', colnames(entry))
colnames(exit) <- gsub('_ex', '', colnames(exit))

# also remove all dashes and underscores in column names
colnames(entry) <- gsub('_', '', colnames(entry))
colnames(exit) <- gsub('_', '', colnames(exit))
colnames(entry) <- gsub('-', '', colnames(entry))
colnames(exit) <- gsub('-', '', colnames(exit))
colnames(ipip) <- gsub('_', '', colnames(ipip))
colnames(pers) <- gsub('_', '', colnames(pers))
colnames(ipip) <- gsub('-', '', colnames(ipip))
colnames(pers) <- gsub('-', '', colnames(pers))

# add rows of entry and exit
d <- rbind(entry, exit)

# keep only rows if they show up twice, so we have a within-measurement of each 
# individual at time t_0 and t_1
d <- ddply(d, .(ID), function(x){
    if(nrow(x)==1){
        return(NULL)
    }
    else{
        return(x)
    }
})

table(d$t)
## 
##  0  1 
## 13 13

We have \(13\) individuals with entry and exit surveys completely filled out. Let’s also add the ipip and pers datasets, and then use only individuals with complete data.

d <- inner_join(d, pers, by = "ID")
d <- inner_join(d, ipip, by = "ID")

After joining these two datasets we have 12 subjects in the dataset and the total sample size is \(n=\) 24.

3 Model designs

3.1 Assumptions concerning the data-generation process

Next, in order to design an appropriate model, one of the things we need to make assumptions about is the underlying data-generative model, i.e., what model generated the type of data we have at our disposal. In short, there are four candidates: Cumulative, Continuation ratio, Stopping ratio, and Adjacent category. Above we opted for a default Cumulative. Let’s instead use LOO to create identical models and then use the same data for all models to compare them from an information theoretical perspective.

bf1 <- bf(mvbind(MAASQ116,
                 MAASQ216,
                 MAASQ316,
                 MAASQ416,
                 MAASQ516,
                 MAASQ616,
                 MAASQ716,
                 MAASQ816,
                 MAASQ916,
                 MAASQ1016,
                 MAASQ1116,
                 MAASQ1216,
                 MAASQ1316,
                 MAASQ1416,
                 MAASQ1516) ~ 
            1 )

m0 <- brm(bf1,
    family = cumulative,
    data = d,
    silent = TRUE,
    refresh = 0
    )
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 17.5 seconds.
## Chain 3 finished in 17.5 seconds.
## Chain 1 finished in 17.8 seconds.
## Chain 4 finished in 18.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 17.7 seconds.
## Total execution time: 18.5 seconds.
mcr <- brm(bf1,
    family = cratio,
    data = d,
    silent = TRUE,
    refresh = 0
    )
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 4.1 seconds.
## Chain 1 finished in 4.1 seconds.
## Chain 3 finished in 5.3 seconds.
## Chain 4 finished in 5.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 4.7 seconds.
## Total execution time: 5.6 seconds.
msr <- brm(bf1,
    family = sratio,
    data = d,
    silent = TRUE,
    refresh = 0
    )
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 5.2 seconds.
## Chain 4 finished in 6.2 seconds.
## Chain 1 finished in 6.3 seconds.
## Chain 3 finished in 6.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 6.1 seconds.
## Total execution time: 6.8 seconds.
mac <- brm(bf1,
    family = acat,
    data = d,
    silent = TRUE,
    refresh = 0
    )
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 8.5 seconds.
## Chain 2 finished in 8.5 seconds.
## Chain 4 finished in 8.6 seconds.
## Chain 3 finished in 9.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 8.6 seconds.
## Total execution time: 9.3 seconds.
m0 <- add_criterion(m0, criterion = "loo")
msr <- add_criterion(msr, criterion = "loo")
mcr <- add_criterion(mcr, criterion = "loo")
mac <- add_criterion(mac, criterion = "loo")

loo_compare(m0, mcr, msr, mac)
##     elpd_diff se_diff
## m0   0.0       0.0   
## mcr -2.7       3.1   
## msr -2.8       3.1   
## mac -3.1       1.6

Evidently, since \(\Delta\)SE is fairly large, in comparison to the relative difference in expected log pointwise predictive density (\(\Delta\)elpd), we can safely assume that there’s no significant difference between the likelihoods and, thus, opt for the standard approach, i.e., the Cumulative distribution, in this case represented by \(\mathcal{M}_0\).

3.2 weekly data

3.2.1 Prior and posterior predictive checks

We want to model the weekly instrument and how it varies over time \(t\). In this case, \(t\) will be modeled using a Gaussian Process and then each individual in the dataset will be modeled using a varying intercept.

weekly <- left_join(weekly, pers, by = "ID")

# make sure ID is a factor
weekly$ID <- as.factor(weekly$ID)

# Need to format data so that each person's measurements goes 
# from 1 ... n, according to order, then we use those numbers as temporal var
weekly <- weekly[order(weekly$ID),]
weekly$Seq <- sequence(tabulate(weekly$ID))
# Pick first question. We could model this as multivariate but it won't give us 
# anything extra
p <- get_prior(WQ1_1_6 ~ gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
               data = weekly,
               family = cumulative())

p$prior[1] <- "normal(0,2)"
p$prior[8] <- "normal(0,5)"
p$prior[13] <- ""
p$prior[14] <- "inv_gamma(1.6, 0.1)"
p$prior[15] <- "weibull(2,1)"

WQ1gp <- brm(WQ1_1_6 ~ gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
             data = weekly,
             family = cumulative(), 
             prior = p, 
             sample_prior = "only",
             silent = TRUE, 
             refresh = 0)
## Running MCMC with 4 chains, at most 8 in parallel...
## 
## Chain 1 finished in 0.2 seconds.
## Chain 3 finished in 0.3 seconds.
## Chain 4 finished in 0.3 seconds.
## Chain 2 finished in 0.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 0.7 seconds.
pp_check(WQ1gp, type="bars", nsamples = 100)

The priors are fairly uniform on the outcome space. Let’s sample with data now.

WQ1gp <- brm(WQ1_1_6 ~ gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
             data = weekly,
             family = cumulative(), 
             prior = p, 
             control = list(adapt_delta = 0.9999, max_treedepth = 13),
             silent = TRUE, 
             refresh = 0)
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 90.4 seconds.
## Chain 1 finished in 91.5 seconds.
## Chain 3 finished in 95.1 seconds.
## Chain 4 finished in 98.7 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 93.9 seconds.
## Total execution time: 99.0 seconds.
WQ2gp <- brm(WQ2_1_6 ~ gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
             data = weekly,
             family = cumulative(), 
             prior = p, 
             control = list(adapt_delta = 0.99, max_treedepth = 12),
             silent = TRUE, 
             refresh = 0)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 71.1 seconds.
## Chain 3 finished in 88.0 seconds.
## Chain 2 finished in 104.7 seconds.
## Chain 4 finished in 106.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 92.6 seconds.
## Total execution time: 106.8 seconds.
WQ3gp <- brm(WQ3_1_6 ~ gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
             data = weekly,
             family = cumulative(), 
             prior = p, 
             control = list(adapt_delta = 0.99),
             silent = TRUE, 
             refresh = 0)
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 47.7 seconds.
## Chain 1 finished in 50.9 seconds.
## Chain 2 finished in 51.2 seconds.
## Chain 4 finished in 53.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 50.7 seconds.
## Total execution time: 53.3 seconds.
WQ4gp <- brm(WQ4_1_6 ~ gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
             data = weekly,
             family = cumulative(), 
             prior = p, 
             control = list(adapt_delta = 0.999),
             silent = TRUE, 
             refresh = 0)
## Running MCMC with 4 chains, at most 8 in parallel...
## 
## Chain 4 finished in 67.8 seconds.
## Chain 3 finished in 74.0 seconds.
## Chain 2 finished in 85.3 seconds.
## Chain 1 finished in 87.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 78.5 seconds.
## Total execution time: 89.0 seconds.
WQ5gp <- brm(WQ5_1_6 ~ gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID),
             data = weekly,
             family = cumulative(), 
             prior = p, 
             control = list(adapt_delta = 0.99, max_treedepth = 12),
             silent = TRUE, 
             refresh = 0)
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 126.8 seconds.
## Chain 4 finished in 130.7 seconds.
## Chain 1 finished in 132.6 seconds.
## Chain 2 finished in 133.6 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 130.9 seconds.
## Total execution time: 134.0 seconds.

3.2.2 Diagnostics

models <- list(WQ1gp, WQ2gp, WQ3gp, WQ4gp, WQ5gp)

# Check divergences, tree depth, energy for each model
for(m in models) {
  rstan::check_hmc_diagnostics(eval(m)$fit)
}
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
# Check rhat and neff
for (m in models) {
  if (max(bayesplot::rhat(eval(m)), na.rm = T) >= 1.01)
    print("Warning: Rhat >= 1.01")
  
  if (min(bayesplot::neff_ratio(eval(m)), na.rm = T) <= 0.2)
    print("Warning: ESS <= 0.2")
}

3.2.3 Parameter estimates

In short, very little happening over time according to the weekly instrument.

3.3 daily data

3.3.1 Prior and posterior predictive checks

Let’s redo this but for the daily dataset.

daily <- left_join(daily, pers, by = "ID")
daily$ID <- as.factor(daily$ID)
daily$date_s <- as.numeric(daily$date)
daily$Living14 <- as.numeric(daily$Living14)
daily <- daily[order(weekly$ID),]
daily$Seq <- sequence(tabulate(daily$ID))

p <- get_prior(resp ~ 1 + gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
               data = daily, family = cumulative)
p$prior[1] <- "normal(0,2)"
p$prior[6] <- "normal(0,5)"
p$prior[17] <- "inv_gamma(1.6, 0.1)"
p$prior[18] <- "weibull(2,1)"

m_daily <- brm(resp ~ 1 + gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
               data = daily, 
               family = cumulative,
               prior = p,
               sample_prior = "only",
               silent = TRUE, 
               refresh = 0)#, threads = threading(4))
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 0.3 seconds.
## Chain 3 finished in 0.3 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 4 finished in 0.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.3 seconds.
## Total execution time: 0.9 seconds.
pp_check(m_daily, type = "bars", nsamples = 100)

m_daily <- brm(resp ~ 1 + gp(Seq) + Age + Gender01 + Occupation + Living14 + (1 | ID), 
               data = daily, 
               family = cumulative,
               prior = p,
               control = list(adapt_delta=0.99),
               silent = TRUE, 
               refresh = 0)
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 67.4 seconds.
## Chain 4 finished in 68.3 seconds.
## Chain 3 finished in 96.6 seconds.
## Chain 1 finished in 103.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 83.8 seconds.
## Total execution time: 103.3 seconds.
pp_check(m_daily, type = "bars", nsamples = 100)

3.3.2 Diagnostics

rstan::check_hmc_diagnostics(m_daily$fit)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
# check rhat and ESS
if(max(bayesplot::rhat(eval(m_daily)), na.rm=T) >= 1.01) 
  print("Warning: Rhat >=1.01")

if(min(bayesplot::neff_ratio(eval(m_daily)), na.rm=T) <= 0.2) 
  print("Warning: ESS <=0.2")

3.3.3 Parameter estimates

Let’s plot an interesting effect, which albeit is not significant has a tendency that’s interesting.

plot(conditional_effects(m_daily), plot = FALSE)[[2]] +
  xlim(1,2.5) +
  scale_x_continuous(name=NULL, 
                     breaks = c(1,2), 
                     label = c("Man/\n transman","Woman/\n transwoman"), 
                     expand = c(0.08, 0)) +
  scale_y_continuous(name=NULL, breaks = c(6,7)) + 
  theme_tufte()

Next, look at the trend.

plot(conditional_effects(m_daily), plot = FALSE)[[5]] +
  xlab("Day") + ylab("") +
  theme_tufte()

One things clear, since the uncertainty explodes at the end (the gray area grows a lot) it makes it very hard to make any inferences. There might be a difference here, but we can’t possibly see it since the uncertainty has so large.

3.4 Multivariate models with temporal variable

We want a model that uses personal data (e.g., age) as predictors. The idea is to later see if any of the predictors have predictive capacity for our six different survey instruments.

First, design a model with predictors from personal data and compare with our null model (\(\mathcal{M}_0\)), together with our indicator \(t\) for time, and a varying intercept ID that varies depending on subject.

bf1 <- bf(mvbind(MAASQ116,
                 MAASQ216,
                 MAASQ316,
                 MAASQ416,
                 MAASQ516,
                 MAASQ616,
                 MAASQ716,
                 MAASQ816,
                 MAASQ916,
                 MAASQ1016,
                 MAASQ1116,
                 MAASQ1216,
                 MAASQ1316,
                 MAASQ1416,
                 MAASQ1516) ~ 
            1 + Age + Gender01 + Occupation + Living14 + t + (1|c|ID))

m_pers <- brm(bf1,
    family = cumulative,
    data = d,
    silent = TRUE,
    refresh = 0)#,
## Running MCMC with 4 parallel chains...
## 
## Chain 2 finished in 175.6 seconds.
## Chain 1 finished in 186.2 seconds.
## Chain 4 finished in 194.9 seconds.
## Chain 3 finished in 198.7 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 188.9 seconds.
## Total execution time: 199.1 seconds.
    #threads = threading(4))

m_pers <- add_criterion(m_pers, criterion = "loo")

# compare the new model with our null model, it should hopefully be better...
(l <- loo_compare(m0, m_pers))
##        elpd_diff se_diff
## m_pers   0.0       0.0  
## m0     -51.3      24.5

Calculating the credible interval (\(z_{95\%}=1.96\)) we can see that it doesn’t cross zero, i.e., adding these predictors makes the relative out of sample prediction capability jump up even though we add seven population-level predictors and one group-level predictor (Living14 contains three parameters to estimate).

l[2,1] + c(-1,1) * 1.96 * l[2,2]
## [1] -87.327351   2.811668

3.4.1 Prior and posterior predictive checks

Let’s now sample all models with the same predictors we used above (visual prior and posterior predictive checks were done for each model).

################################################################################
#
# MAAS
#
################################################################################
bf1 <- bf(mvbind(MAASQ116,
                 MAASQ216,
                 MAASQ316,
                 MAASQ416,
                 MAASQ516,
                 MAASQ616,
                 MAASQ716,
                 MAASQ816,
                 MAASQ916,
                 MAASQ1016,
                 MAASQ1116,
                 MAASQ1216,
                 MAASQ1316,
                 MAASQ1416,
                 MAASQ1516) ~ 
            1 + Age + Gender01 + Occupation + Living14 + t + (1 |c| ID))

# Prior predictive checks have been conducted and the probability mass 
# is distributed evenly on the outcome space
p <- get_prior(bf1, data=d, family=cumulative) %>%
  mutate(prior = ifelse(class == "b", "normal(0,3)", prior)) %>%
  mutate(prior = ifelse(class == "cor", "lkj(2)", prior)) %>%
  mutate(prior = ifelse(class == "Intercept", "normal(0,5)", prior)) %>%
  mutate(prior = ifelse(class == "sd", "weibull(2,1)", prior))
  
         
m_maas <- brm(bf1,
    family = cumulative,
    data = d,
    prior = p,
    #sample_prior = "only", # use if you only want to sample from prior
    silent = TRUE,
    refresh = 0)#,
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 31.2 seconds.
## Chain 1 finished in 31.3 seconds.
## Chain 2 finished in 31.8 seconds.
## Chain 3 finished in 32.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 31.7 seconds.
## Total execution time: 32.6 seconds.
    #threads = threading(4))

# use the below to do prior and posterior predictive checks
# change resp to the response variable you want to check 
# pp_check(m_maas, resp = "MAASQ116", type = "bars, nsamples = 250) 
################################################################################
#
# SPANE
#
################################################################################

bf1 <- bf(mvbind(SPANEQ115,
                 SPANEQ215,
                 SPANEQ315,
                 SPANEQ415,
                 SPANEQ515,
                 SPANEQ615,
                 SPANEQ715,
                 SPANEQ815,
                 SPANEQ915,
                 SPANEQ1015,
                 SPANEQ1115,
                 SPANEQ1215) ~ 
            1 + Age + Gender01 + Occupation + Living14 + t + (1 |c| ID))

p <- get_prior(bf1, data=d, family=cumulative) %>%
  mutate(prior = ifelse(class == "b", "normal(0,3)", prior)) %>%
  mutate(prior = ifelse(class == "cor", "lkj(2)", prior)) %>%
  mutate(prior = ifelse(class == "Intercept", "normal(0,5)", prior)) %>%
  mutate(prior = ifelse(class == "sd", "weibull(2,1)", prior))

m_spane <- brm(bf1,
    family = cumulative,
    data = d,
    prior = p,
    # sample_prior = "only", # use if you only want to sample from prior
    silent = TRUE,
    refresh = 0)#,
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 21.5 seconds.
## Chain 4 finished in 21.7 seconds.
## Chain 2 finished in 21.9 seconds.
## Chain 1 finished in 22.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 21.8 seconds.
## Total execution time: 22.4 seconds.
    #threads = threading(4))

# use the below to do prior and posterior predictive checks
# change resp to the response variable you want to check 
# pp_check(m_spane, resp = "SPANEQ115", type = "bars", nsamples = 250) 
################################################################################
#
# PWB
#
################################################################################

bf1 <- bf(mvbind(PWBQ117,
                 PWBQ217,
                 PWBQ317,
                 PWBQ417,
                 PWBQ517,
                 PWBQ617,
                 PWBQ717,
                 PWBQ817) ~ 
            1 + Age + Gender01 + Occupation + Living14 + t + (1 |c| ID))

p <- get_prior(bf1, data=d, family=cumulative) %>%
  mutate(prior = ifelse(class == "b", "normal(0,3)", prior)) %>%
  mutate(prior = ifelse(class == "cor", "lkj(2)", prior)) %>%
  mutate(prior = ifelse(class == "Intercept", "normal(0,5)", prior)) %>%
  mutate(prior = ifelse(class == "sd", "weibull(2,1)", prior))

m_pwb <- brm(bf1,
    family = cumulative,
    data = d,
    prior = p,
    # sample_prior = "only", # use if you only want to sample from prior
    silent = TRUE,
    refresh = 0)#,
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 18.0 seconds.
## Chain 1 finished in 18.4 seconds.
## Chain 2 finished in 18.4 seconds.
## Chain 3 finished in 19.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 18.5 seconds.
## Total execution time: 19.5 seconds.
    #threads = threading(4))

# pp_check(m_pwb, resp = "PWBQ117", type = "bars", nsamples = 250) 
################################################################################
#
# PST
#
################################################################################
# Here we move to a bernoulli since we have 0/1 outcome
bf1 <- bf(mvbind(PSTQ101,
                 PSTQ201,
                 PSTQ301,
                 PSTQ401,
                 PSTQ501,
                 PSTQ601,
                 PSTQ701,
                 PSTQ801,
                 PSTQ901,
                 PSTQ1001,
                 PSTQ1101,
                 PSTQ1201,
                 PSTQ1301,
                 PSTQ1401,
                 PSTQ1501,
                 PSTQ1601,
                 PSTQ1701,
                 PSTQ1801,
                 PSTQ1901,
                 PSTQ2001,
                 PSTQ2101,
                 PSTQ2201) ~ 
            1 + Age + Gender01 + Occupation + Living14 + t + (1 |c| ID))

p <- get_prior(bf1, data=d, family=bernoulli) %>%
  mutate(prior = ifelse(class == "b", "normal(0,3)", prior)) %>%
  mutate(prior = ifelse(class == "cor", "lkj(2)", prior)) %>%
  mutate(prior = ifelse(class == "Intercept", "normal(0,5)", prior)) %>%
  mutate(prior = ifelse(class == "sd", "weibull(2,1)", prior))

m_pst <- brm(bf1,
    family = bernoulli,
    data = d,
    prior = p,
    # sample_prior = "only",
    silent = TRUE,
    refresh = 0
    )
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 11.0 seconds.
## Chain 3 finished in 11.1 seconds.
## Chain 2 finished in 11.2 seconds.
## Chain 4 finished in 11.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 11.2 seconds.
## Total execution time: 11.7 seconds.
# pp_check(m_pst, resp = "PSTQ101", type = "bars", nsamples = 250) 
################################################################################
#
# SE - two models (Cumulative and Bernoulli), both multivariate 
#
################################################################################
# NOTE: the first two questions can be analyzed as Bernoulli 
bf1 <- bf(mvbind(#SEQ114, # this and the next outcome only has two categories
                 #SEQ214,
                 SEQ314,
                 SEQ414,
                 SEQ514,
                 SEQ614,
                 SEQ714,
                 SEQ814,
                 SEQ914,
                 SEQ1014) ~
            1 + Age + Gender01 + Occupation + Living14 + t + (1 |c| ID))

p <- get_prior(bf1, data=d, family=cumulative) %>%
  mutate(prior = ifelse(class == "b", "normal(0,3)", prior)) %>%
  mutate(prior = ifelse(class == "cor", "lkj(2)", prior)) %>%
  mutate(prior = ifelse(class == "Intercept", "normal(0,5)", prior)) %>%
  mutate(prior = ifelse(class == "sd", "weibull(2,1)", prior))

m_se <- brm(bf1,
    family = cumulative,
    data = d,
    prior = p,
    # sample_prior = "only",
    silent = TRUE,
    refresh = 0)#,
## Running MCMC with 4 parallel chains...
## 
## Chain 4 finished in 8.0 seconds.
## Chain 1 finished in 8.4 seconds.
## Chain 2 finished in 9.4 seconds.
## Chain 3 finished in 9.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 8.8 seconds.
## Total execution time: 9.8 seconds.
    #threads = threading(4))

# pp_check(m_seq, resp = "SEQ314", type = "bars", nsamples = 250) 

# next take the two questions that only had answers on two levels
bf1 <- bf(mvbind(SEQ114,
                 SEQ214) ~
            1 + Age + Gender01 + Occupation + Living14 + t + (1 |c| ID))

p <- get_prior(bf1, data=d, family=bernoulli) %>%
  mutate(prior = ifelse(class == "b", "normal(0,3)", prior)) %>%
  mutate(prior = ifelse(class == "cor", "lkj(2)", prior)) %>%
  mutate(prior = ifelse(class == "Intercept", "normal(0,5)", prior)) %>%
  mutate(prior = ifelse(class == "sd", "weibull(2,1)", prior))

m_se2 <- brm(bf1,
    family = bernoulli,
    data = d,
    prior = p,
    # sample_prior = "only",
    silent = TRUE,
    refresh = 0
    )
## Running MCMC with 4 parallel chains...
## 
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 finished in 0.4 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.7 seconds.
# pp_check(m_seq2, resp = "SEQ114", type = "bars", nsamples = 250) 
################################################################################
#
# PPHQ
#
################################################################################


bf1 <- bf(mvbind(PPHQ115,
                 PPHQ215,
                 PPHQ315,
                 PPHQ415,
                 PPHQ515,
                 PPHQ615,
                 PPHQ715,
                 PPRQ1110,
                 PPRQ2110,
                 PPRQ3110,
                 PPOQ117
                 ) ~ 
            1 + Age + Gender01 + Occupation + Living14 + t + (1 |c| ID))

p <- get_prior(bf1, data=d, family=cumulative) %>%
  mutate(prior = ifelse(class == "b", "normal(0,3)", prior)) %>%
  mutate(prior = ifelse(class == "cor", "lkj(2)", prior)) %>%
  mutate(prior = ifelse(class == "Intercept", "normal(0,5)", prior)) %>%
  mutate(prior = ifelse(class == "sd", "weibull(2,1)", prior))

m_pph <- brm(bf1,
    family = cumulative,
    data = d,
    prior = p,
    # sample_prior = "only",
    silent = TRUE,
    refresh = 0,
    # threads = threading(4)
    )
## Running MCMC with 4 parallel chains...
## 
## Chain 3 finished in 18.5 seconds.
## Chain 2 finished in 18.8 seconds.
## Chain 4 finished in 19.2 seconds.
## Chain 1 finished in 19.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 19.0 seconds.
## Total execution time: 19.8 seconds.
# pp_check(m_pphq, resp = "PPHQ115", type = "bars", nsamples = 250)

So we have now sampled seven models. For each instrument we have one multivariate model (\(>1\) outcome), where we also estimate the correlation between the outcome given the subject’s ID.

Most of the models are using, as we’ve already discussed, the Cumulative family. However, in a few cases we use the Bernoulli family. First, we use it for outcomes that were binary (i.e., Yes/No outcomes). Second, we use it for two questions in the SE instrument (\(\mathcal{M}_{\text{SE}2}\)), which only had empirical outcomes on two levels.

Before we put any trust in these models we need to check a number of diagnostics.

3.4.2 Diagnostics

First, we check general HMC diagnostics. Second, we check that \(\widehat{R} \leq 1.01\) and that the effective sample size (ESS) is \(\geq 0.2\).

models <- list(m_maas, m_spane, m_pwb, m_pst, m_se, m_se2, m_pph)

# Check divergences, tree depth, energy for each model
for(m in models) {
  rstan::check_hmc_diagnostics(eval(m)$fit)
}
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
# Check rhat and neff
for (m in models) {
  if (max(bayesplot::rhat(eval(m)), na.rm = T) >= 1.01)
    # diag in corr matrix always NA
    print("Warning: Rhat >= 1.01")
  
  if (min(bayesplot::neff_ratio(eval(m)), na.rm = T) <= 0.2)
    print("Warning: ESS <= 0.2")
}

Hamiltonian Monte Carlo diagnostics indicate all is well (one can also check traceplots for all models, i.e., using plot(model)).

3.4.3 Parameter estimates

First, let’s summarize what we have so far:

  1. Weekly trends modeled with Gaussian Processes (results in Sect. 3.2).
  2. Daily trends modeled with Gaussian Processes (results in Sect. 3.3). We’ve also shown conditional effects of gender since it is the main driver for differences we’re seeing.
  3. Seven multivariate models with temporal variable to analyze entry/exit. Additionally, we have an eight model for analyzing two questions (Q1 and Q2) in the SE instrument (using a Bernoulli likelihood).

We can now investigate which, if any, parameters are significant on the 95%-level for the third item in the list, i.e., for each model, and each question, plot all parameters. Any parameters that are significant could warrant further analysis.

For each model we first analyze the temporal variable \(t\). Then we analyze the rest of the \(\beta\) parameters. Finally, where appropriate we look at the conditional effects. All densities we plot are limited to 95% probability mass, while the 50% probability mass is shaded around a vertical line, which is the median.

In short, any density that does not cross zero (the vertical line) is potentially a significant effect.

Each of these models used a sample size of \(n=\) 24 for sampling.

3.4.3.1 Mindful attention awareness scale (MAAS)

The model estimated 1355 parameters.

For Questions \(1\)–\(3\), \(8\), and \(14\) the parameter \(t\) is negative. In short, the respondents answered with lower responses at \(t_1\) than at \(t_0\) for these questions. If we look at each unique parameter, in each question, can we see which, if any, predictors drove these results?

Much variance makes the estimates uncertain and, ultimately, no parameter is significant (not even Q7, even though it looks that way).

Two notable parameters, first for Q1 (woman/transwoman answers significantly higher than man/transman), while in Q6 man/transman answers significantly lower.

Occupation status (student/non-student) does not have any effect. Let’s, finally, look at living conditions before we turn our attention to the next instrument.

For Q15, if living with a partner, the respondent will respond significantly higher. The same for Q13 if the respondent lives with their family.

3.4.3.2 Scale of positive and negative experience (SPANE)

3.4.3.3 Personal well-being index (PWB)

For Q4 the respondents gave higher responses at \(t_1\). Let’s look at the parameters for that question.

No parameters influencing that question to a larger extent. In short, respondents, no matter their age, etc., responded higher at \(t_1\) for Q4.

3.4.3.4 Positive thinking (PST)

Nothing here.

3.4.3.5 Self efficacy (SE)

Nothing here.

Nothing here.

3.4.3.6 Perceived productivity (PP)

Questions \(2\), \(6\), and \(11\), show a significant difference when moving from \(t_0\) to \(t_1\) (lower responses at \(t_1\)). Let’s examine the parameters for these three questions.

The predictor age has a negative effect (the higher the age the lower the response). How does it look like as a conditional effect? A conditional effect keeps all predictors fixed except for the one we’re interested in, i.e., age.

conditional_effects(m_pph, effects="Age", resp = "PPHQ215")

So age seems to be the main driver and this is how it looks like. Zero on the \(x\)-axis is in this case the median age, i.e., approximately \(36\) years.

Gender has a negative effect (male/transmale answers lower on the response), however it is actually not significant. The only thing we can say is that a number of parameters makes respondents, on average, respond with lower values at \(t_1\).

Living with a partner makes respondents answer with higher response in Q11, however it is not significant since it does cross zero. The only thing we can say is that a number of parameters makes respondents, on average, respond with lower values at \(t_1\).

4 Computational environment

print(sessionInfo(), locale=FALSE)
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] latex2exp_0.4.0   ggthemes_4.2.4    here_1.0.1        data.table_1.14.0
##  [5] openxlsx_4.2.3    patchwork_1.1.1   bayesplot_1.8.0   forcats_0.5.1    
##  [9] stringr_1.4.0     dplyr_1.0.4       purrr_0.3.4       readr_1.4.0      
## [13] tidyr_1.1.2       tibble_3.1.0      ggplot2_3.3.3     tidyverse_1.3.0  
## [17] plyr_1.8.6        brms_2.14.11      Rcpp_1.0.6       
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1            backports_1.2.1         igraph_1.2.6           
##   [4] splines_4.0.4           crosstalk_1.1.1         TH.data_1.0-10         
##   [7] rstantools_2.1.1        inline_0.3.17           digest_0.6.27          
##  [10] htmltools_0.5.1.1       rsconnect_0.8.16        fansi_0.4.2            
##  [13] magrittr_2.0.1          modelr_0.1.8            RcppParallel_5.0.3     
##  [16] matrixStats_0.58.0      xts_0.12.1              sandwich_3.0-0         
##  [19] prettyunits_1.1.1       colorspace_2.0-0        rvest_0.3.6            
##  [22] haven_2.3.1             xfun_0.21               callr_3.5.1            
##  [25] crayon_1.4.1            jsonlite_1.7.2          lme4_1.1-26            
##  [28] survival_3.2-7          zoo_1.8-8               glue_1.4.2             
##  [31] gtable_0.3.0            emmeans_1.5.4           V8_3.4.0               
##  [34] pkgbuild_1.2.0          rstan_2.26.0.9000       abind_1.4-5            
##  [37] scales_1.1.1            mvtnorm_1.1-1           DBI_1.1.1              
##  [40] miniUI_0.1.1.1          xtable_1.8-4            stats4_4.0.4           
##  [43] StanHeaders_2.26.0.9000 DT_0.17                 htmlwidgets_1.5.3      
##  [46] httr_1.4.2              threejs_0.3.3           ellipsis_0.3.1         
##  [49] farver_2.1.0            pkgconfig_2.0.3         loo_2.4.1              
##  [52] sass_0.3.1              dbplyr_2.1.0            utf8_1.1.4             
##  [55] labeling_0.4.2          tidyselect_1.1.0        rlang_0.4.10           
##  [58] reshape2_1.4.4          later_1.1.0.1           munsell_0.5.0          
##  [61] cellranger_1.1.0        tools_4.0.4             cli_2.3.1              
##  [64] generics_0.1.0          broom_0.7.5             ggridges_0.5.3         
##  [67] evaluate_0.14           fastmap_1.1.0           yaml_2.2.1             
##  [70] processx_3.4.5          knitr_1.31              fs_1.5.0               
##  [73] zip_2.1.1               nlme_3.1-152            mime_0.10              
##  [76] projpred_2.0.2          xml2_1.3.2              compiler_4.0.4         
##  [79] shinythemes_1.2.0       rstudioapi_0.13         curl_4.3               
##  [82] gamm4_0.2-6             reprex_1.0.0            statmod_1.4.35         
##  [85] bslib_0.2.4             stringi_1.5.3           highr_0.8              
##  [88] ps_1.6.0                Brobdingnag_1.2-6       lattice_0.20-41        
##  [91] Matrix_1.3-2            nloptr_1.2.2.2          markdown_1.1           
##  [94] shinyjs_2.0.0           vctrs_0.3.6             pillar_1.5.0           
##  [97] lifecycle_1.0.0         jquerylib_0.1.3         bridgesampling_1.0-0   
## [100] estimability_1.3        httpuv_1.5.5            R6_2.5.0               
## [103] bookdown_0.21           promises_1.2.0.1        gridExtra_2.3          
## [106] codetools_0.2-18        boot_1.3-27             colourpicker_1.1.0     
## [109] MASS_7.3-53.1           gtools_3.8.2            assertthat_0.2.1       
## [112] rprojroot_2.0.2         withr_2.4.1             shinystan_2.5.0        
## [115] multcomp_1.4-16         mgcv_1.8-34             parallel_4.0.4         
## [118] hms_1.0.0               grid_4.0.4              coda_0.19-4            
## [121] minqa_1.2.4             cmdstanr_0.3.0.9000     rmarkdown_2.7          
## [124] shiny_1.6.0             lubridate_1.7.10        base64enc_0.1-3        
## [127] dygraphs_1.1.1.6

5 References

  1. https://ppc.sas.upenn.edu/resources/questionnaires-researchers/mindful-attention-awareness-scale↩︎

  2. https://www.midss.org/content/scale-positive-and-negative-experience-spane-0↩︎

  3. http://www.acqol.com.au/instruments↩︎

  4. https://www.psytoolkit.org/survey-library/pts.html↩︎

  5. https://sparqtools.org/mobility-measure/new-general-self-efficacy-scale/↩︎

  6. https://www.hcp.med.harvard.edu/hpq/info.php↩︎

  7. https://ipip.ori.org/MiniIPIPKey.htm↩︎

  8. https://www.psykiatri-regionh.dk/who-5/who-5-questionnaires/Pages/default.aspx↩︎